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Abstract 


The Helmholtz equation is solved within a three-dimensional rectangular duct with a 
sound source at the duct entrance plane, local admittance conditions on the side walls, 
and a new, nonlocal radiation boundary condition at the duct exit plane. The formulation 
employs a truncation of an infinite matrix, the generalized modal admittance tensor , that 
represents the transformation of the modal pressure coefficients to the modal axial velocity 
coefficients. This condition accurately models the acoustic admittance (a relationship be- 
tween the velocity and pressure of an acoustic wave) at an arbitrary computational boundary 
plane. In particular, the proper physical admittance may correspond to a nonreflecting radi- 
ation condition on an arbitrarily located boundary plane that is used to limit the size of the 
computational domain. A linear system of equations is constructed with second-order central 
differences for the Helmholtz operator and second-order backward differences for both local 
admittance conditions and the gradient term in the nonlocal radiation boundary condition. 
The resulting matrix equation is large, sparse, and non-Hermitian. The size and structure 
of the matrix makes direct solution techniques impractical; as a result, a nonstation ary it- 
erative technique is used for its solution. The theory behind the nonstationary technique 
is reviewed, and numerical results are presented for radiation from both a point source and 
a planar acoustic source in both soft-walled and hard-walled ducts. The solutions with the 
nonlocal boundary conditions are invariant to the location of the computational boundary, 
and the same nonlocal conditions are valid for all solutions. The nonlocal conditions thus 
provide a means of minimizing the size of three-dimensional computational domains. 




1 Introduction 


Computational aeroacoustics requires the resolution of relative short waves in large compu- 
tational volumes. In the case of aircraft engines, for example, the frequencies of interest may 
require a grid spacing of roughly one centimeter and a corresponding computational element 
volume of one cubic centimeter. Since modern engines have volumes of several cubic meters, 
a computational aeroacoustic simulation must involve millions of equations. 

Acoustics problems are often posed in an Infinite domain, but a computational domain 
must be finite. The computational domain is made finite by the introduction of a computa- 
tional boundary condition on it’s arbitrarily-located surface within the infinite domain. Most 
computational boundary conditions for these surfaces are local differential operators which 
presume some knowledge of the wave field in the neighborhood of the boundary surface, [1], 
[ 2 ]. 

If the surface is far from the sources of sound, this is a good approximation because 
theoretical acoustics gives a good estimate of the direction of wave propagation. But moving 
the bounding surface away from the sources makes the computional problem large. Reducing 
the size of the computational problem by moving the boundary surface close to the sources 
introduces errors because the complex wave fields near the source are incompatible with the 
assumptions used in the construction of the boundary conditions. 

We are developing computational aeroacoustic boundary conditions which are valid for 
all wave fields, irrespective of their complexity. These boundary conditions may be applied 
close to the acoustic sources, so that the size of the computational volumes may be kept 
to a minimum. These conditions represent the superposition of all possible linear fields in 
the domain exterior to the computational volume. Since linearity is the only assumption 
about the exterior field, the computational volume may be reduced to the limit of this 
linear assumption. These boundary conditions are naturally more complicated because they 
contain more information. They are implemented as matrix or tensor operators and are called 
nonlocal boundary conditions because they represent the influence of a wave propagation 
between separate points on the boudary surface. 

The first formulation and demonstration of these nonlocal boundary conditions was given 
in reference [3] for the case of wave propagation in a two-dimensional duct. There, it was 
shown that the computational space could be reduced by a a factor of ten with no effect on 
the solution. A comparison of the nonlocal boundary condition with some local boundary 
conditions in reference [4] showed that the nonlocal conditions were accurate for all source 
fields, whereas the local conditions were sometimes accurate and sometimes inaccurate. 

The purpose of this paper is to extend the formulation of reference [3] to the case of 
a three-dimensional sound field. Again, we use a rectangular duct so that the boundary 
surface is a simple plane. Because we must consider interactions between all pairs of points 
on this boundary surface, the boundary condition is a tensor operator. 

The following section describes the problem in detail. Then the tensor boundary condition 
is derived for the case of a semi-infinite duct. This is followed by a discrete formulation and 
solution process for the Helmholtz (wave) equation with the nonlocal boundary condition. 
Computational results are given for several point sources oscillating at both low and high 
frequencies, and these results are compared to computations from convergent series solutions. 
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Outgoing Wave Section 



Figure 1. Semi- infinite duct geometry. 


2 Problem Description 

We wish to demonstrate the utility of a nonlocal boundary condition on a simple three- 
dimensional acoustics problem with a known solution. To this end, we consider the semi- 
infinite rectangular duct shown in figure 1. We want to find the numerical solution of the 
Helmholtz equation wfith a monopole source: 

V 2 p + k 2 p — s (1) 

with local admittance boundary conditions on the entrance boundary plane, at x = 0, and 
on the side walls. These local admittance boundary conditions have the form 


V n = 


r z0P + Vs 
pc 


( 2 ) 


w r here p and c. are the ambient density and sound speed, respectively. This problem arises 
in duct acoustics and has a series solution. Physically, p is the complex acoustic pressure in 
the duct, and k = uijc , where u is the frequency of the disturbance. That is, 


p{x,t) = p(x)e ,u " 


( 3 ) 


where i is the unit complex number 1. The physical side conditions are specified by 
the local dimensionless admittances /?, and a nonlocal admittance radiation condition will 
be utilized on the surface x = X. This nonlocal boundary condition is derived from the 
following two constraints: 

1. The general solution in the “outgoing wave” section (x > I), including the boundary 
(x = X), can be given as a series of radiating acoustic modes, or eigenfunctions, which 
satisfy the side w r all boundary conditions. 

2. The solution and its derivative with respect to x are continuous across the computa- 
tional boundary plane (x = X), as required by the axial component of the momentum 
equation and by the continuity equation at the boundary plane. 

For the problem considered here, these conditions lead directly to a natural nonlocal bound- 
ary condition. 
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Real duct problems are neither semi-infinite or likely to be limited to a finite set of modes, 
given an arbitrary source, so that the boundary condition will be only approximate. However, 
the intention here is to formulate a discrete boundary condition suitable for coupling with an 
internal discretization. In finite discretizations, only a finite number of modes are possible, 
and the discrete boundary condition will couple naturally with the interior discretization. 
Note that in an experimental setting, data at a radiation boundary are naturally described 
in terms of exterior modes and we show that this boundary condition is also useful in a 
numerical setting. 

The discrete radiation boundary condition is a tensor equation of the form 

= T zftjkj'k'Pj'k' T &sjk (4) 

pc 

where pyw is the complex pressure at grid point (j\ k f ), Ujk is the axial velocity at grid point 
(j, &), and repeated indices are summed (Einstein summation notation) over th e finite range 
of integers j and k. The index 5 refers to the source position, while all other indices are related 
to the numerical grid. For example, they correspond to samples of the pressure and axial 
velocity at grid points j/ ; and on planes where x is constant. Boundary condition (4) is a 
generalization of the two-dimensional nonlocal boundary condition developed in reference [3]. 
Reference [4] showed that this nonlocal boundary condition is more accurate than several 
other local boundary conditions in the two-dimensional case. 

The nonlocal boundary condition is computationally more expensive than local boundary 
conditions; however, section 4 provides evidence that the nonlocality is necessary to specify 
a proper physical interface between the computational domain and the exterior domain at 
an arbitrary boundary plane in a duct. Local boundary conditions involve some assumptions 
about the waves and their direction of propagation at the boundary. They do not permit 
general wave fields in the exterior domain. 

The three-dimensional discretization of the Helmholtz equation with the nonlocal bound- 
ary conditons results in a large, sparse non-Hermitian matrix equation for the pressures in 
the computational domain. The matrix dimension for the higher frequency problems was of 
order 10 . Because of the structure and large size of this matrix, attempts at solutions with 
direct methods were unsuccessful in the high-frequency cases. Solutions were attempted with 
two iterative, nonstationary Krylov subspace techniques: the QMR and TFQMR methods 
outlined in section 4. These iterative methods succeeded after the matrix was precondi- 
tioned. These techniques proved to be quite robust and are used to get the results shown in 
section 5. 
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3 Tensor Admittance Boundary Condition 

Assume that at the boundary plane x = L the velocity data are linearly related to the 
pressure data by an admittance tensor 


1 a 

u{L,yj,z k ) = —0jkj'k'P(L : yj',z k >) (5) 

pc 

where u is the axial velocity and p is the pressure. The tensor Pjkj'k' will be referred to as 
the nodal admittance tensor. In general, the tensor fd depends on the general solution in the 
exterior domain. We may, for example, specify that the exterior solution is composed only 
of modes that describe right w r ard moving waves in the duct in figure 1. In many practical 
instances, known analytic forms of the exterior solution exist. Radiating modes from an 
infinite duct will be used here for analysis. 

Note that calculation of the modes is not an explicit part of the solution scheme even 

though they are used to construct the boundar}^ condition tensor /3j k ji k r. This construction 
is done only once in advance of any calculations in the computational domain, and the 
boundary condition data are read or linked into the computational solution scheme. The 
following section gives a detailed formulation of the nonlocal admittance tensor for the infinite 
three-dimensional duct radiation condition. 


3.1 Infinite Three-Dimensional Duct Radiation Condition 

Modal solutions to rectangular-duct acoustic equations are 

?£»(*»»»*) = pC 2 <f>mn(y,z)e ±tkxmn{X ~ L) ( 6 ) 

The function <j> mn {y, z) is the mode function and the time factor e~ fWf is understood. Further, 
the sign of the imaginary part of the complex axial wavenumber k xmn is positive. The mode 
function is defined by 

tpmniy ♦ c) == Xm (0 

Xm(y) = a m e+ tk ^y + b m e- tk ^y (8) 

0n(*) = c n e +tkznZ -f d n e~' k * nZ (9) 

krmn = yjk* - + k* m ) ( 10 ) 

where eigenvalues k ym and k zn are computed with the methods in appendix A and the modal 
constants are computed using the methods of appendix B. In the simplest case of a hard- 
walled duct, k ym — mir/H and k zn = mr/W. In soft-walled ducts, these eigenvalues are 
complex variables. 

The pressure and velocity in a field of progressive waves in the infinite duct for (x > L ) 
is 


p(x,y,z) = pc 2 ' P mn4>mn{y,z)e +tkxmn{X (11) 

771 ,71 = 0 
OO 

u{x,y,z) = c ^ U+ n ({> mn (y,z)e +tkimn(x - L) (12) 

m,n=0 


where V+ n and U+ n are arbitrary constants that represent the modal pressure and velocity 
amplitudes. This solution is the general wave field solution in the exterior domain. The 
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axial momentum equation gives an alternate formula for the velocity 


\(x,y,z) = c 


oo / r 
ft'X 


m, n=0 


k 


KntmniV, z)e + ' k ^ x ~ L) 


( 13 ) 


Consistency of the two formulas requires that the velocity and pressure coefficients be related 
by 


p+ = f t ^ u+ 
mn \ 1 77171 
\ fcxmn / 


The constants of proportionality are called the modal impedances Z mn \ 

k \ 


Zm n — 


k x 


(14) 


(15) 


With these modal impedances, the proportionality of pressure and velocity is expressed as 

T^rnn = ( 16 ) 

The above impedance formula is a special adaptation of a general formula that relates the 
coefficients of any source-free wave field. 


^mn — ^ ^ ^mnm t n r 

m',n r = 0 


(17) 


Equation (17) is a general modal impedance boundary condition that can represent an} 7 linear 
homogeneous boundary condition for the duct. It can also represent radiation effects at the 
end of a finite duct. The modal boundary condition that represents radiation in an infinite 
duct is recovered from the general formula with the definition 


&mm t &nn r 


(18) 


where 8 mm ' is the Kronecker delta. The general modal formulas for the pressure and velocity 
at the plane x = L are, then, 


P{k, yiZ'j — pc ) ’ 4 > mn{yiZ') ^ ) Z mnm i n i ld m i n ‘ 

m,n=0 


u 


{L,y,z) = c 4>mn{y- l z)U„ 

m,n=0 


(19) 

( 20 ) 


3.2 Nodal Impedance Boundary Condition 

Modal solutions to duct acoustic problems are effective when the duct is uniform or has a 
finite number of uniform sections joined together [5]. In regions where the duct is highly 
nonuniform, general numerical methods are needed. Examples of nonuniform regions are 
transitions around centerbodies and sections that contain turbomachinery. Numerical meth- 
ods applied to these regions need boundary conditions expressed in terms of the data struc- 
tures of the method — typically, the values of the dependent variables at grid points, or nodes , 
of the numerical method. Assuming that the dependent variables are pressure and velocity 
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and that the acoustic field exterior to the computational domain is defined by linear equa- 
tions, then a linear relation must exist between these variables on the boundary surface, 
which is part of both the interior, or computational, domain and the exterior linear domain. 
With a two-dimensional surface grid of points ( j/y , z*), the most general linear homogeneous 
relation between pressure and velocity over the grid of nodes is called the nodal impedance 

Zjkj'k 1 - 

P]k = PC Zjkj'k'Uj'k' (21) 

Exterior sources may be included by adding an inhomogeneous term to the above equation: 

Pjk — PC ^ . Z jkj' k'^j'k 1 T Psjk (22) 

i'J '= 1 

The unspecified summation ranges in the nodal impedance formulas are finite. We assume a 
rectangular computation grid with 1 < j < J, 1 < k < K for simplicity. The points yj,z k 
are taken to be equally spaced with intervals Ay,Az. These assumptions can be relaxed; 
however, this step adds nothing to the general implementaton here. 

A relationship clearly exists between the modal impedance and the nodal impedance. 
The relation is made explicit by utilizing the nodal velocities to solve for the modal velocity 
coefficients and placing the resulting solution in an expression for the nodal pressures. The 
solution, in its simplest form, requires an equal number of modes and nodes in each coordinate 
dimension. 

One- dimensional data structures are needed for the two-dimensional arrays. At a fixed 
axial location (x = T), a vector of nodal velocities is defined by 

u jk = u(L,y } ,z k ) (23) 

{ujk} = {{fife},} (24) 

where the second script, varies through its range before the first script j is incremented. 
The vector on the left side of equation (23) is a column vector of column vectors. A similar 
structure is adopted for the nodal pressures: 

Pit - p(L,Vi,z k) (25) 

{M = {{«},} < 26 ) 

The nodal velocities are related to the modal velocity coefficients through the finite matrix 
approximation of the infinite sum in equation (20): 

{u jk } = c{(j> jkrnn ]{U mn } (27) 

tfrjkmn — 4 > mn{yj'i %k) (^) 

In the above matrix, the first two scripts, jk are row dependent and the second two scripts 
are column dependent. 

The modal velocity coefficients are given by a matrix inversion 

{ u mn } = c" 1 [d> ]k mn ] -1 { u jk } (29) 

Now, the nodal pressure is the finite matrix approximation to the infinite sum of equa- 
tion (19): 

{ Pjk } = P^ [^jlmn] [-^mnm'n'] { ld m ' n ' } 

= pC [Zmnm’n'] [4 > j'k'm'n'} { Hj'k' } (30) 
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Consequently, the nodal impedance is given by the following similarity transformation of the 
modal impedance : 

[Zjkj'k'] = (31) 

This result gives the desired nodal impedance relation between pressure and velocity: 

{ Pjk } = pc [Zjkj'k 1 ] { Uj'k 1 } (32) 


3.3 Nodal Admittance Boundary Condition 

Other boundary conditions are readily formed with the nodal impedances. For example, a 
Helmholtz solver that utilizes only pressure as the dependent variable would require a relation 
between the pressure and the pressure gradient on the bounding surface. This relation must 
be nonlocal in order to be completely general. 

In the present case of a duct for which the computational domain is bounded by the 
plane x = L, the axial momentum equation gives the axial velocity as 

= (33) 

The pressure gradient and pressure are then related by 

{(EbJ = '*(/WH;y*'} (34) 

where the nodal admittance matrix is the inverse of the nodal impedance matrix: 


[0jkj'k>] = [Zjkj'k 1 ] 


-1 


(35) 


This nodal admittance matrix is clearly a similarity transformation of a modal admittance 
matrix. Elements of this modal admittance matrix for an infinite duct would be the inverse 
of the elements in equation (18). Let the modal admittance be defined as the inverse of the 
modal impedance as shown: 

[Bmnm'n' ] ~ [Zmnm'n’] 1 (36) 

The radiation condition for an infinite uniform duct is 


B 


mnm’n' — 



(37) 


The nodal admittance is a similarity transformation of the modal admittance : 




-1 


(38) 


3.4 Tensor Product Transformations 

The matrices that define the boundary condition transformations become large for high- 
frequency computations. If we have a duct lateral dimension of 2 3 = 8 wavelengths and the 
same number of points per wavelength in the computation, then each dimension requires 
2 6 = 64 grid points. The impedance and admittance matrices would have a dimension of 
2 12 = 4096. Direct inversion should be avoided for these large matrices. Fortunately, the 
matrices in the boundary conditions are tensor products so that the computations are done 
separately for each dimension with the simplifications of the tensor product identities. These 
identities are given here to show the explicit computation procedures. 
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The similarity matrix can be written as 



— [X;m[0Arn]] 

(39) 

[Xjm] 

= [Xm(]/;)] 

(40) 

[0Jbn] 

= [li>n(zk)] 

(41) 


The similarity matrix [<f>jkmn] is formed with blocks of the matrix [xpkn] scaled by elements 
of the matrix [Xm(Pj)]* Tensor product notation for this form is 

= [Xjm] ® (42) 

In this notation, the second matrix supplies the block that is scaled by an element of the 
first matrix. With this convention, the inverse of the tensor product is the tensor product 
of the inverses, taken in the same order: 

[(ijimn] = [Xjm] ® [t^fcn] 

The admittance and impedance matrices transform a vector into another vector. Let x, y 
be ( N x 1) matrices and A be an (N x N) matrix. The transformation y = Ax then takes 
N multiplications per row and N additions per row, assuming y is initially zero. Thus, a 
single transformation has 2N 2 operations. The admittance transformation is 

{{«}} = M ® [0] l B ) [XT 1 ® [^] _1 { {?} } (44) 

Each subvector is an (A r x 1) matrix, and each sub-block matrix is an (N x A r ) matrix. The 
transformation can proceed from right to left. The first operation is to multiply [0] -1 times 
each of the N subvectors {p}. This step requires 2N 3 operations, and returns N subvectors. 
The next operation is to take linear combinations of these subvectors with the elements 
of each row of the matrix [x]” 1 as coefficients. This step also requires 2N 3 operations and 
produces a vector of N 2 elements. Premultiplication of this vector by the diagonal matrix [B] 
requires A r 2 operations. Each of the next two operations takes 2A^ 3 operations, so that the 
total number of operations is 8A r 3 + N 2 for each transformation. By contrast, the full form 
of the matrix transformation would require 2(A r2 ) 2 = 2A M operations. The tensor product 
form is superior to the expanded form for computation with iterative techniques. The tensor 
product is also superior for storage. Each matrix has N 2 elements, so that the total storage 
is 5A r2 words. The expanded form of the transformation would use N 4 words. 


(43) 
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4 The Discrete System 

The discrete system represents the Helmholtz equation and all boundary conditions in dis- 
crete form. The Laplacian operator in three dimensions is represented by a seven-point 
stencil as follows: 

1 

V 2 p,jfc = a,iij>k'Pi+i' < j+j>,k+k' (45) 

i'j’k '=- 1 

where the stencil data depend on the algorithm that is implemented. Data for the 

second-order- accurate central-difference scheme are 


a 0,0,0 — 

- 2 l 


1 

a ±l,0,0 = 

Aa: 2 


1 

a 0,±l,0 = 

Ay 2 


1 

a 0,0,±l = 

Az 2 


1 

A^2 


+ 


1 L 
Ay 2 + 



(46) 

(47) 

(48) 

(49) 


where Ax, Ay, and A z are the grid spacings in the respective dimensions. 

Local admittance boundary conditions on the “entrance plane” and sidewalls are all taken 
in the form of equation (2). The momentum equations for the x, y, and 2 directions are 


dp , — 

— — |- ikpcu 
ox 

dp , — 
— — h ikpcv 
dy 


dp 

dz 


+ ikpcw 


= 0 
= 0 
= 0 


(50) 

(51) 

(52) 


If we combine equations with the local admittance boundary conditions on the entrance and 
side walls, then we obtain 


dp i,j,k 
dx 


ikpc0oj t kPij,k 


dpixk 

dy 

dpi,j,k 

dy 

dpi,j,l 

dz 


— lkpC0_ff/2,i,kPi,l,k 

+ lkpC0+n/2,i,kPi,J,k 

— lkpC0-W/2,ijPi,j,\ 


dpij.K 

dz 


T lkpC0+W/2,i,jPi,j,K 


ikpcugjk , 1 < j < J, 1 < k < K 

(53) 

1 <i< I, 

1 <k<K 

(54) 

1 < i < I, 

1 < k < I< 

(55) 

1 <i<I, 

1 <J<J 

(56) 

1 < i < I, 

1 <J<J 

(57) 


A discrete nonlocal admittance boundary condition, as described in the previous section, is 
used on the exit plane x = xj — L. 


dfirpyc 

dx 




(58) 
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Figure 2. Sparse matrix from the linear system. Nonzero elements are indicated by dark 
marks. The block in the lower right corner is generated by the nonlocal boundary condition. 


The sum in equation (58) is over the entire y-z grid at the boundary. Forward or back- 
ward difference formulas were utilized at the computational domain surfaces to avoid the 
introduction of “ghost points.” 

The three-dimensional array of unknowns is arranged as a column matrix, where 
the indices i, j , and k vary in reverse order with their position; that is, k varies first, then j 
varies, and then i varies last. 


{ Ptjk } = P(X«, Vj, Zk) = { { { Pk }j }, } (59) 

{u, jk } = U(Xi,yj,Zk) = { { { «Jt }j } ( 60 ) 

Major partitions of this column of unknowns are based on the axial index i. Each partition 
is a column matrix that contains the unknown pressures on a single cross-sectional plane. 
Subpartitions of these columns are based on the y-dimension index j. These subpartitions 
define column matrices for which only the z-dimension index k varies. 

The system matrix is the square coefficient matrix of the above column of unknowns. 
Its nonzero elements are the coefficients from the difference operator (45)— (49) , the local 
boundary conditions (53)-(57), and the nonlocal admittance boundary condition (58). The 
system matrix is very sparse (most of its elements are zero) because the Helmholtz equation 
is a local operator; however, a significant block of (J x K) nonzero elements exists due to 
the nonlocal boundary condition. Figure 2 is a graphic illustration of the system matrix 
structure. This system matrix is clearly sparse, except for the large data block in the lower 
right corner where the coefficients for the nonlocal boundary condition are stored. 

Direct solutions to the above system matrix were attempted and found to be impractical, 
at least for workstation calculations. Banded matrix solvers could not be used because they 
fill the inner null bands of the sparse matrix with nonzero elements. With matrix dimensions 
0[1O 5 ], the matrix bandwidth would be 0 [10 3 ], and a banded solver would create O [10 8 ] 
elements. Consequently, the iterative QMR and TFQMR techniques described below were 
selected to solve the discrete system matrix equation. These methods are known as Krylov 
subspace techniques. Further details on these techniques can be found in references [6], 
[7], and [8]. In accordance with the recommendation of Freund, [9], the Krylov subspace 
techniques were applied to the system matrix in its complex form. Otherwise, when a 
complex system of dimension N is written as a real system of dimension 2A r , then a matrix 
may be produced with a distribution of eigenvalues that includes the origin or values very 
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near the origin. A matrix with such values is singular or nearly singular and can present 
difficulties for the QMR and TFQMR Krylov subspace techniques. 

4.1 Preconditioning Concepts 

In some cases, the system matrix equation could not be solved with the iterative techniques 
without the application of preconditioning. This preconditioning was necessary for conver- 
gence of the iteration (not just for acceleration of the rate of convergence). In other cases 
that could be solved without preconditioning, the preconditioning accelerated the rate of 
convergence. Thus, preconditioning was always beneficial. 

The notion of preconditioning is a simple. Instead of solving the original system Ax = b, 
an equivalent system A'x' = b' is solved, where A' = PAQ, x r = Q~*x, and b' = Pb. The 
equivalent system must be constructed such that the matrix A' has more desirable properties 
than the original matrix A. An example of a more desirable property is an improved matrix 
condition number, which makes the equivalent system amenable to solution by iterative 
techniques. 

If an invertable matrix can be chosen in factored form M = M\M 2 , then the elements of 
the equivalent preconditioned system A'x' = b' are of the form 

A' = M 1 ~ l AM 2 -\ b' = Mr 1 6, x' = M 2 x 

The corresponding Krylov subspace technique generates in a translated subspace, 

x' n € x' 0 + )C n (r' 0 , A') (62) 

where JC n (r' 0 , A') = span {r(,, A'r' Q , A'^r' 0: ■ ■ • , A ,n-1 r(,}. In terms of the original system, 
Krylov iterates x n and residuals r n = b — Ax n are connected by 

x n = M 2 ~ l x' n € xo + IC n (M~ l ro, M~ l A) and r„ = Mir' n (63) 

For minimum residual techniques, the right preconditioning matrix M 2 = I is the identity. 
This identity matrix is generally preferred so that the preconditioned residual is the same as 
the ordinary residual. 

The more common preconditioners are reviewed in reference [10]. The choice of precon- 
ditioners used here was influenced more by their availability than by a detailed evaluation 
of their suitability for the problem at hand. The ones used here were members of a general 
class of preconditioners that are implemented within the software package qmrpack [11]. 



4.2 Methods of Preconditioning 

Two preconditioning methods were found that worked effectively. These methods were the 
Symmetric Successive Overrelaxation (SSOR) method and the Incomplete Gaussian Elimi- 
nation (ILU) method. These methods are described more fully below. 

4.2.1 SSOR Method 

The SSOR method is an interative technique for solving a linear system Ax = b. The matrix 
A is written as D — E — F, where D , E , and F are the diagonal, strict lower triangular, and 
strict upper triangular portions of A , respectively. For the classic Successive Overrelaxation 
(SOR) technique, a parameter u, which is the overrelaxation parameter , is introduced and 
Ax = b is rewritten as 

( D — uE)x — [(1 — <jo)D + ujF] x = cab (64) 
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An iterative method is produced by letting the leftmost x be the new iterate id m+l) and the 
rightmost x be the old iterate 

(D - uE)u (rn+x] = [(1 - u)D + uF] u {m) + ub (65) 

or 

ytm+i) _ p _ [(j _ + u (m) + ( D _ uE y l ub ( 66 ) 

The inverse is only slightly more expensive to compute than a matrix multiplication because 
the matrix is lower triangular. 

The SSOR technique is a modification of the SOR technique in which the vector is 
modified from the top to the bottom and then from the bottom back up to the top to 
complete one iteration. The SSOR technique has the form 

u (m+i/2) _ (£>_ w £)-i[(i - u )D + u>F]ul m) + {D-uE)- l u>b 
u (m+i) = (£ _ uF)- 1 [(1 -u)D+uE\u< m+1/2) + (D-uE)- l ub (67) 

Written in one step, SSOR is an iteration u( m+1 * = MvS m ^ + /, where 

M = I — - 1 )(D - uF)~ 1 D(D - u>E)~ x A (68) 

and / is a fixed vector that may be worked out from (67). The idea in any iterative technique 
is to minimize M in some way (usually by minimizing the spectral radius). Ideally, M = 
0 so that the inversion can be accomplished in one step. If M is “small,” then (D — 
uF)~ l D(D — t oE)~ l is a good approximate inverse of A (modulo a constant multiplier); 
moreover, because of the triangular form of its factors, it is easily computed. This matrix 
M is the preconditioner in the SSOR technique. 

4.2.2 ILU Method 

The ILU method has variants that are designated by ILU(n), where n is an integer. The 
idea that underlies the basic ILU(O) preconditing is taken from Gaussian elimination. First, 
find an upper triangular matrix U and lower triangular L such that A — LU. The ILU 
preconditioning contains a remainder term A = LU + R and imposes a constraint on the 
structure of L -f U (i.e., L -f U must have nonzero elements within the envelope in which A 
has nonzero elements). The remainder term results because fill-in is not allowed at elements 
of L and U that were already zero. With this constraint, the underlying matrix-vector 
multiplication with the preconditioned matrix is no more expensive than when applied to 
the original matrix A. The ILU(O) may be extended by allowing additional inward fill in L 
and U, which leads to the ILU(p) algorithms, where p is the maximum inward fill. Further 
details of the ILU(O) algorithm and its variants may be found in [10]. 

Numerical experiments with both preconditioning methods showed that they were effec- 
tive in accelerating the rates of convergence of the iterative solutions. The ILU method was 
found to be more reliable for this problem; however, so that it w r as adopted to prepare the 
results shown in the following section. 
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5 Results 


This section contains the results of numerical solution of the discrete system outlined in 
section 4. The computationally less difficult results are presented first and the more difficult 
results are presented later. The higher frequencies and more complicated sources present 
increased difficulties. Two types of sources are specified at the inlet boundary that represent 
opposite extremes: a plane wave source and a point source. The source conditions are 
specified as 

^ f - = v, jk 1 < j < J, 1 <k<K 

If v B jk — C, where C is a complex constant, the source is a plane wave source. If = 



Figure 3. Approximate representation of the delta function by an inverted cone. Support of 
the cone extends over five grid points. 

C S(y — yo, z — zq), the source is a point source. Series solutions for both source boundary 
conditions are derived in appendix C. The delta- function source is not representable by a 
single continuous function or on a discrete grid. In the following numerical experiments, the 
truncated modal series was used to analytically represent the delta function. Similarly, an 
inverted cone of unit volume, such as the one in figure 3, was used to numerically represent 
the delta function. 



Figure 4. Duct computational domain. The shaded region indicates the horizontal plane on 
which pressure data are shown in subsequent plots. 

The plotted results to be shown later are the complex pressure on a horizontal plane that 
cuts through the center of the duct as shown in figure 4. Some plots show the result of the 
application of the boundary condition at different locations in the duct. The purpose here 
is to show that the solution is fairly insensitive to the point of application of the boundary 
condition. In each computation, a uniform grid with equal grid spacing in each coordinate 
direction was used. The grid spacing was choosen to give a resolution of 10 points per 
wavelength. 
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5.1 Hard- Walled Duct with Plane Source 

In a hard- walled duct, the admittance of each of the four side walls is identically zero. Plots 
of the real part of acoustic pressure are shown in Figure 5 for a planar source that oscillates 
at a frequency of 500 Hz. As expected, the acoustic pressure is independent of the spanwise 
direction and no evidence of any reflections from the radiation boundary is observed. The 
interior solution is in excellent agreement with the analytical solution on this scale. Real 
and imaginary parts of the acoustic pressure differed by a phase shift of 90°, so only the real 
part of the acoustic pressure is plotted here. 
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(a) Series. (b) Numerical. 

Figure 5. Real part of the pressure field due to a 500 Hz plane-w T ave source in a hard-walled 
duct. Results computed with a (17 x 17 x 17) grid in a 1.0 m cubic domain. 

Figure 6 shows the results obtained for the real part of the acoustic pressure, where the 
frequency of the planar source has been increased to 1,000 Hz. Here, the number of grid 
points per meter was doubled, to accommodate the shorter wavelength, and the radiation 
boundary condition was applied at 0.5 m down the duct (preserving the total number of grd 
points). Results are consistent with those obtained at 500 Hz. Again, the imaginary part of 
the acoustic pressure was not plotted because it was identical to the real part, except for a 
phase shift of 90°. 



5.2 Hard- Walled Duct with Point Source 

To demonstrate that the solution is independent of the radiation boundary condition for 
all sources, a point source was used to generate the acoustic wave pattern in the hard wall 
duct. Point source calculations are a severe test of both the numerical method and radiation 
boundary condition. Classic modal analysis shows that eigenfunction components of the 
acoustic fields generated by point sources are of two types: modes that decay exponentially 
in space (cut-off modes) and modes that are purely periodic in space (cut-on modes). The 
real and imaginary parts of the acoustic pressure calculated for a point source that oscillates 
at 1,000 Hz are shown in figures 7 and 8 respectively. The radiation boundary condition 
was applied 1.0 m from the source. Several cut-on modes generate a complicated pressure 
pattern away from the source. Near the radiation boundary, the acoustic pressure field, 
computed from the numerical method, is identical to the modal series results. However, in 
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(a) Series. 


(b) Numerical. 


Figure 6. Real part of the pressure field due to a 1000 Hz plane- wave source in a hard-walled 
duct. Results computed with a (17 x 33 x 33) grid in a section with axial dimension of 0.5 m. 


the vicinity of the point source, the imaginary part of the acoustic pressure computed from 
the modal series and the numerical method show a discrepancy. The discrepancy between 
the numerical method and modal series in the vicinity of the point source was expected. 
Near the source plane, a large number of duct modes and, hence, grid points are needed to 
represent the acoustic point source. Thus, the discrepancy near the source is believed to be 
caused by the coarseness of the grid near the source. 

To further investigate the discrepancy in the computed acoustic field near the point source 
at 1,000 Hz, numerical results were obtained when the radiation boundary was moved to 
0.5 m from the source, These results are shown in figures 9 and 10. 

Note that the movement of the radiation boundary closer to the source had little or 
no effect on the acoustic pressure field computed by the numerical method. This result 
lends further credence to the hypothesis that the discrepancy near the source is due to 
the coarseness, of the grid (because the grid spacing was not changed when the radiation 
boundary was moved closer to the source). The basic solution process would be expected to 
work on an uneven grid; however, this case was not tested. 

For point sources that oscillate at frequencies above 1000 Hz, the combined effect of many 
cut-on and cut-off modes generates a pattern that indicates multidirectional wavefronts that 
w T ould not be captured with simpler local-plane-wave type of boundary conditions. Numerical 
results are plotted in figures 11 and 12 for a point source that oscillates at 1,500 Hz. The 
radiation boundary condition was applied at 1.0 m from the source. Again, a discrepancy 
can be observed in the imaginary part of the pressure in the vicinity of the point source. 
When the radiation boundary is placed at 0.5 m from the source (figures 13 and 14, this 
discrepancy is still present. 
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x, m x, m 

(a) Series. (b) Numerical. 


Figure 7. Real part of the pressure field due to a 1000 Hz point source in a hard- walled duct. 
Nonlocal boundary condition applied at 1.0 m. Results computed with a (33 x 33 x 33) grid. 
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(a) Series. (b) Numerical. 


Figure 8. Imaginary part of the pressure field due to a 1000 Hz point source in a hard-walled 
duct. Nonlocal boundary condition applied at 1.0 m . Results computed with a (33 x 33 x 33) 
grid. 
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(a) Series. (b) Numerical. 

Figure 9. Real part of the pressure field due to a 1000 Hz point source in a hard-walled duct. 
Nonlocal boundary condition applied at 0.5 m. Results computed with a (17 x 33 x 33) grid. 



(a) Series. 


(b) Numerical. 


Figure 10. Imaginary part of the pressure field due to a 1000 Hz point source in a hard-walled 
duct. Nonlocal boundary condition applied at 0.5 m. Results computed with a (17 x 33 x 33) 
grid. 
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(a) Series. (b) Numerical. 

Figure 11. Real part of the pressure field due to a 1500 Hz point source in a hard- walled 
duct. Nonlocal boundary condition applied at 1.0 m. Results computed with a (49 x 49 x 49) 
grid. 
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(a) Series. (b) Numerical. 


Figure 12. Imaginary part of the pressure field due to a 1500 Hz point source in a hard-walled 
duct. Nonlocal boundary condition applied at 1,0 m. Results computed with a (49 x49 x49) 
grid. 
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Figure 13. Real part of the pressure field due to a 1500 Hz point source in a hard- walled 
duct. Nonlocal boundary condition applied at 0.5 m. Results computed with a (25 x 49 x 49) 
grid. 
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(a) Series. (b) Numerical. 

Figure 14. Imaginary part of the pressure field due to a 1500 Hz point source in a hard-walled 
duct. Nonlocal boundary condition applied at 0.5 m. Results computed with a (25 x 49 x 49) 
grid. 
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5.3 Soft-Walled Duct with Point Source 

The radiation condition has also been tested with acoustically treated, or soft, walls. In 
contrast to hard-walled ducts, acoustically treated ducts have surfaces that are not perfectly 
reflecting. Results for the soft-walled duct case were computed with the admittance j3 = 
.5 — .5 1 in equation (2). Acoustically treated walls damp out acoustic disturbances and reduce 
the number of contributing waves compared with the hard- walled case. To limit the number 
of graphs, results for the soft-walled duct are only presented for the more-troublesome point 
source. 

Figures 15 through 18 show the results of a point source that oscillates at 1,000 Hz in 
the soft-walled duct. Figures 15 and 16 show results for the radiation boundary located 
1.0 m from the source; figures 17 and 18 show results for the radiation boundary located 
at 0.5 m from the point source. On this scale, the numerical solution is not affected by the 
distance of the radiation boundary from the point source. The discrepancy that occurred in 
the imaginary part of the acoustic pressure for hard-walled ducts is not observed in the lined 
duct. This difference is believed be caused by the fact that the lining has a damping effect on 
the acoustic waves and also reduces the number of contributing waves in comparison with to 
those in the rigid-w r alled duct. The agreement between the results for the numerical solution 
and modal series is excellent. Reflections from the nonlocal radiation boundary condition 
are not evident at this frequency. 

Figures 19 through 22 show results of a 1500 Hz point source that oscillates in the lined 
duct. Significant damping of the acoustic w T aves is clearly observed, in comparison with the 
corresponding results for the hard-w r alled duct. The conclusions drawn from this set of plots 
are consistent with those for 1,000 Hz. 
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(a) Series. (b) Numerical. 


Figure 15. Real part of the pressure field due to a 1000 Hz point source in a soft-walled duct. 
Nonlocal boundary condition applied at 1.0 m. Results computed with a (33 x 33 x 33) grid. 
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(a) Series. (b) Numerical. 


Figure 16. Imaginary part of the pressure field due to a 1000 Hz point source in a soft-walled 
duct. Nonlocal boundary condition applied at 1.0 m. Results computed with a (33 x 33 x 33) 
grid. 
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Figure 17. Real part of the pressure field due to a 1000 Hz point source in a soft- walled duct. 
Nonlocal boundary condition applied at 0.5 m. Results computed with a (17 x 33 x 33) grid. 
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(a) Series. (b) Numerical. 

Figure 18. Imaginary part of the pressure field due to a 1000 Hz point source in a soft-walled 
duct. Nonlocal boundary condition applied at 0.5 m. Results computed with a (17 x 33 x 33) 
grid. 
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Figure 19. Real part of the pressure field due 
Nonlocal boundary condition applied at 1.0 
grid. 
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(b) Numerical. 


a 1500 Hz point source in a soft-walled duct. 
. Results computed with a (49 x 49 x 49) 
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(b) Numerical. 


Figure 20. Imaginary part of the pressure field due to a 1500 Hz point source in a soft-walled 
duct. Nonlocal boundary condition applied at 1.0 m. Results computed with a (49 x 49 x 49) 
grid. 
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(a) Series. 


(b) Numerical. 


Figure 21. Real part of the pressure field due to a 1500 Hz point source in a soft-walled duct. 
Nonlocal boundary condition applied at 0.5 m. Results computed with a (25 x 49 x 49) grid. 
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(a) Series. (b) Numerical. 

Figure 22. Imaginary part of the pressure field due to a 1500 Hz point source in a soft-walled 
duct. Nonlocal boundary condition applied at 0.5 m. Results computed with a (25 x 49 x49) 
grid. 
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6 Concluding Remarks 

A nonlocal boundary condition has been formulated for the Helmholtz equation in three- 
dimensional Cartesian coordinates. The boundary condition represents the nonreflection 
condition on a cross-sectional plane within a semi-infinite rectangular duct. This boundary 
condition is completely independent of the wave field and transmits waves from any source 
distribution without perceptible reflection. The boundary plane can be moved close to the 
source without affecting the solution for the wave field. 

In order to demonstrate the boundary condition, a three-dimensional system of finite- 
difference equations was written to model acoustic waves generated by arbitrary sources in 
a semi-infinite duct with the Helmholtz equation. Source conditions were prescribed at the 
entrance plane, and conventional local boundary conditions were used on the duct sidewalls. 
The nonlocal boundary condition was applied on the exit or end plane of the duct, which 
defines the extent of the computational domain. On this end plane, which can be placed at 
an arbitrary position, the acoustic field must be the same as the field within a semi-infinite 
duct with an arbitrary source distribution. 

The resulting system of equations was large and sparse. In order to solve this system for 
practical frequencies for which the duct dimensions were many multiples of the wavelength, 
two types of preconditioning were investigated. These two methods were Symmetric Succes- 
sive Over-relaxation (SSOR) preconditioning and Incomplete Gaussian Elimination (ILU) 
preconditioning. The ILU method was found to be more reliable for this problem. System 
solutions were found with ILU preconditioning and Krylov subspace iterative schemes. 

To assess the effectiveness of the nonlocal boundary condition, numerical solutions were 
compared with classical series solutions. The physical model used in the comparison was a 
duct that measured 1 m x 1 m in cross section. The duct was semi-infinite, but the com- 
putational domain was closed at the arbitrarily located right end by the nonlocal boundary 
condition. Comparisons were made for two different source types (i.e., plane waves and point 
sources) at the left end. Comparisons were also made for both hard and absorbing side walls. 
Computations were made for 500, 1,000, and 1,500 Hz at standard atmospheric conditions 
to generate complex sound fields. Within this context, the closing boundary was moved to 
various positions to verify that the predicted sound field was invariant to its location. In 
all cases, the boundary location had no perceptible effect on the predicted sound field. The 
numerically predicted field matched the field produced by the series solution with a high 
degree of accuracy. Because the boundary condition represents the exterior domain general 
solution, it can be used to obtain computional solutions for irregular interior domains when 
classic series solutions are not available. 

These examples demonstrate that nonlocal boundary conditions on the surface of a three- 
dimensional computational domain can accurately represent the effect of a surrounding infi- 
nite domain on the acoustic solutions within the finite computational domain. Because the 
solutions are invariant to the location of the computational boundary, the nonlocal boundary 
conditions provide a means of minimizing the size of the computational domain. 
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A Eigenvalue Computation 


A.l High-Order Eigenvalues 

When the walls are hard, the eigenvalues are multiples of 7r; when the walls have finite 
admittances, a complex characteristic equation must be solved for the eigenvalues. This step 
can easily be the most difficult part of the computation, so the method is developed here in 
considerable detail. We group the admittance $ at the wall with the frequency parameter k 
and the height H by defining To = pcfiokH and tjj = pcfi i- 

If we assume elementary solutions of the form e ±,is,y and apply admittance boundary 
conditions at the lower and upper walls, the following equation is obtained for the eigenvalues 
k y H: 




' i r ££_’ 

1 kyH 

1 -I- JDL - 

1 + kyH 


1 - 


kyH 


1 + 


= 0 


(69) 


(70) 


In the case of hard walls, where r 0 = 77/ = 0, solutions to equation (69) are 

kyH — mir , m — 0,1,2 ,--- 

where S m is used here as a measure of the displacement of an eigenvalue from its nominal 
hard- walled position in the complex plane. 

This elementary case suggests the definition 


(k y H) m = mir + <$ m , m = 0, 1, 2, 


(71) 


A. 1.1 Characteristic Equation for Higher Modes, m > 0 

The characteristic equation for the higher modes can be given in logarithmic form as 

F(S m ) = 6 m 


+ 


5 log 
5 log 


1 + 


1 - 


To 


mir + <5 to J 
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rriK + 6„ 


+ 2 l0g 

4 l0g 


1 + 


1 - 


th 


m-K + 8 m J 
th 

m 7r + 6 m . 


= 0 


(72) 

This form of the characteristic equation shows clearly that the solution 8 m approaches zero 
in the limit where m — » oo. The points at which mi + £ m = ±io,f/ a re singularities of F 
that cannot be solutions. These points should not cause difficulty, except possibly in the 
case where m = 0. The explicit form of equation (72) is important. By computing each 
logarithim separately, we avoid crossing the branch cut, which may introduce an error of 7r 
in the value of the complex log. 


A. 1.2 Newton’s Method for Solution of Higher Modes 

Newton’s method for solution of equation (72) utilizes an iteration for 6 m given by the 
replacement operation 

£ . £ /70\ 

K ” F'(6 m ) (73) 

where the derivative of the characteristic function is 

To 


F\8 m ) = 1-1 


— I- 


TH 


(m7r + <5 m ) 2 — To (mir + 6 m ) 2 - rjf 


(74) 
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This solution method is effective for higher modes (m > 1). One can quickly generate a 
thousand or more eigenvalues on a modern workstation. The work in developing a useful 
basis for a modal computation is limited then to computing the lower eigenvalues. This 
computation is discussed below. 


A. 2 Low-Order Eigenvalues 

A. 2.1 Alternate form of the characteristic equation 

Newton’s method is fast and dependable for the case in which m > 1; however, special care 
is needed for cases in which m = 0 or m = 1. The problem is worse when the admittances 
have large imaginary parts and small real parts, so we attack these cases first. In these 
cases, pure imaginary solutions may exist for the eigenvalue k y H . Accordingly the following 
definitions are introduced to investigate these cases: 

k y H = -ir) (75) 

a = ~t 0 t h (76) 

l t 0 + t h (nn ^ 

b = -i — (7/) 

These changes of variable allow the characteristic equation to be expressed as 

F(tj) = a + r] 2 + 2 brj coth rj = 0 (78) 

To this point, no real change has occurred, because all variables are complex. Now consider 
solutions in which the parameters a b are real variables (i.e., the Teal parts of the complex 
parameters defined above). We can find real solutions £ for these real parameters and add 
the effects of their imaginary parts later. 

A. 2. 2 Solutions for Real Parameters 

Now take the real parts ai b\ of the parameters a = ai + iai and b = b\ + z &2 • For brevity, we 
use the same symbols a b with the understanding that they are restricted to real values in 
this subsection. When the parameters are real, a better form of the characteristic equation 


utilizes the square of the eigenvalue as 



£ = 

v 2 = - ( W 2 

(79) 
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« + { + 2i S(0 

(80) 

o 

Ill 

cosh T] 

(81) 

III 

sinh(77) 

(82) 


The value of the characteristic function at the origin and its behavior when £ is large are 
important for determining if the characteristic function has zeros for positive £. 

F(0) = a + 2b (83) 

F(0 ~ £, £ -> °c (8 4) 

The defined functions C(£) and S(£) have properties that make them useful in determining 
solutions to the characteristic equation. These properties will be listed here before consid- 
ering the numerical procedure for solving for the eigenvalues. 
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The functions are defined by the series for the hyperbolic functions: 

f" 


c(0 = E 


5(0 = E 


e 


n=o (2n + 1)! 

The series definitions permit the following immediate observations: 
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(92) 


The following formulas are useful for evaluating the derivatives of the functions: 


c '(0 = \su) 
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<7(0 - SjQ 
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Finally, the functions have the following asymptotic character: 
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The first derivative of the characteristic function is 
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( 100 ) 


The derivative F'(£) approaches unity from above or below, depending on the sign of b when 
£ is large. The derivative F'( 0) may be positive or negative, depending on the value of 5; 
F'(0) is positive for b > —3/2, zero for b = -3/2, and negative for 6 < -3/2. The second 
derivative of the characteristic function is 


HO = 

b 

[C(0(5 2 (f)-2) + 5(f)l 

(101) 
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The function within brackets is positive, so that the sign of the second derivative depends 
only on the sign of b. If b is negative, then the curvature is positive; if b is positive, the 
curvature is negative. The limit of the function within brackets is 4/45 for £ = 0, which 
gives the value of the derivative at the origin. 

Small solutions |£| << 1 are possible for certain combinations of the parameters a b. 
These combinations are identified by utilizing a two-term series for the characteristic func- 
tion: 

F(£) = F( 0) + F'( OK = 0 (104) 

which gives the following estimate for T ] 2 


t = ~ 


fa + 2b\ 

(rrp) 


(105) 


In order for this estimate to be consistent with the assumed smallness of £, the parameters 
must be related by 

- |F'(0)|e < F(0) < + |F'(0)|e (106) 


1 + f» 


e < F(0) < + 




(107) 


In the special case for which b — > —3/2, the limits above show that a —* —2b. These solutions 
can be positive or negative because they represent the effect of the parameters in moving 
the solutions from the real K y H axis to the imaginary K y H axis in the neighborhood of the 
origin. 

A long list of cases needs to be investigated in determining whether solutions for K y H 
exist near the origin or on the imaginary axis (£ is positive). The cases will be organized here 
by defining ranges of the characteristic function and its slope at the origin. These parameters 
(F(0) and F'(0)) depend on the parameters a b only. 


1. IfF'(0) ^ ~ f- c , then the slope of the characteristic function is positive everyvrhere. If 
F(0) > +c, then no positive solution exists; however, 

( a ) If |F(0)| < e, then a small solution £ = O(e) exists, or else 

(b) else If F(0) < — e, then a positive solution £ > 0 exists, or else 

2. If |F'(0)| < e, then the slope at the origin is near zero. If F(0) > +e 2 , then no positive 
solution exists; however, 

(a) If |F(0)| < e 2 , then a small double solution ^ = O(e) exists, or else 

(b) If F( 0) < — e 2 , then a positive solution exists, or else 

3. If F'(0) < — c, then the slope at the origin is negative and the characteristic function 
has a minimum for some £ > 0. Find the location of the minimum £ m ,„ and evaluate the 
characteristic function to get its minimum value F min . If F mtn > e 2 , then no solution 
exists; however, 

(a) If | F m in | < e 2 , then a double solution exists near £ m i„, or else 

(b) If F m i„ < — e 2 , then the solutions depend on the value of F at the origin. 

i. If F(0) > +e 2 , then two distinct solutions £ < £ m ,- n and £ > £ mt - n exist, or else 

ii. If |F(0)| < e 2 , then a solution near £ = 0 exists and another solution £ > £ mtri 
exists, or else 

iii. If F(0) < — e 2 , then a single solution £ > £ m ; n exists. 
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A. 2. 3 Solutions for Complex Parameters 

Now that we have found all possible solutions for real parameters, the next step is to find 
the effect of the imaginary parts of the parameters on these solutions. Let the parameters 
a b be complex functions of a real parameter 5 , which varies from zero to unity (0 < s < 1). 
The parameters are identical to their real parts when 5 = 0 and take their actual values 
when 5=1. Begin with the solutions for real parts of the parameters for the case in which 
5 = 0 and trace these solutions to the point at which 5 = 1. Note that these solutions can 
include effects of the real parts of the admittances in the parameter a. The tracing follows 
a differential equation for the square of the eigenvalue £($) as a function of the parameter s: 


K,(s)H = y-««) 

(108) 

a = —t 0 th 

(109) 

, ■ To + T H 

6 = ‘ 2 

(110) 

a(s) = 3fc[a] + ?Q[a]s = ai + ia 2 s 

(111) 

b(s) = 3?[6] + i^s[b]s = bi + ib 2 s 

(112) 


Now, £(5), <2(5), and 6(5) are complex functions of the real variable 5. The characteristic 
function depends on s and is complex but is otherwise unchanged from the one defined for 
the real variables: 


m°) 

dms) 

d£ 

&ms) 

oe 

om f) 

ds 


a { s ) + £ + 2^(5) 

( 


1 + b( s) 

M 

s 3 (0 [ 


CTO 

s(0 

(C(()S(0- 


L ) 


mo 

c(0(s 2 (0-0 + s(0' 


2£ 2 


a^s) + 2 b'(s) 


<TO 

s<0 


(113) 

(114) 

(115) 

(116) 


The solutions for which 5 = 0 are 
which is used to find £(5): 


initial conditions for a first-order differential equation, 



(117) 


This equation is valid as long as the initial condition is not a double eigenvalue. When the 
initial condition is a double eigenvalue £ m , the partial derivative in the denominator is 


dF 

d( 


br- 


K - U) 


(118) 


The differential equation for £ is then singular, but a differential equation for the square of 
the displacements of the eigenvalues from the double eigenvalue is not singular. 


d(t-U ) 2 



(119) 
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Given an initial step As, the two eigenvalues near the double eigenvalue are 


( = im ± 



(120) 


The differential equation for the eigenvalue is then used to trace each of these to the final 
value for s = 1. 
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B Acoustic Modes 

The velocity and pressure due to a single acoustic mode are given by 

k 

Admittance boundary conditions on the top and bottom walls are 



r 

CD 

ISO 

1 

Ip J 

. o pc 2 J [ 


k y e tk >' y -V‘* vV 
ke tk * y ke~' k » y 


y ^ ’ — /VyU - - J <3 1 

. Uj 


n -aj{;} = o 
L 1 +AJ ( .) =0 


P. 

We can combine these equations ot obtain a matrix equation for the modal constants: 

( kyh — kh(3 t )e +,k » h / 2 — ( k y h + kh/3 t )e~ ,kvh ^ 2 1 f a 1 _ 

(k y h + kh/3 b )e-' k *V 2 ~{k y h - khp b )e +tk * h ' 2 \ U J “ 


( 121 ) 

( 122 ) 

(123) 

(124) 


B.l Ratio of Coefficients 

Either of the above equations can be used to solve for the ratio a/6 or 6/a. The first gives 


a) {k y h + kh/3 t ) 


(125) 


and the second gives 

/a\ = {kyh - khfol kyh 

\bJ {kyh + khfo) 

The product of these results is unity, which is the eigenvalue equation treated in appendix A. 
Here, however, we want only the ratio, with the assumption that the eigenvalue solution is 
complete. Because some error may occur in the eigenvalue calculation, we take the arithmetic 
average of the two solutions. 





(127) 


B.2 Normalization 


Given the ratio of the coefficients, the magnitude of one coefficient (e.g. ,a) can now be 
chosen such that the average potential energy in the mode is unity. With the two-sided 
transform, the average potential energy is the average of the magnitude of the pressure. 
Consequently, we set 

i r H / 2 7 

(128) 

This condition determines the magnitude of a: 



Way J 


(e ,k » y , 



e ,k * y ) 

,e ,k * y ) 



(129) 
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(130) 


where the inner product of two variables is defined as 

1 f H / 2 

iu ' v) =HLl2 U ‘ Vdy=l 

The inner normalization condition determines the magnitude of a. A choice is still available 
for the phase. The convention will be used that the phase of the pressure is zero at y = 0: 

argp(0) = 0 — > arg(a + b) = 0 (131) 

This choice is achieved by setting 

arg a = — arg [1 + ( bja)\ = 0 (132) 


B.3 Inner Products 

The inner products needed are 


The sine function is defined with a limiting case 


sin (feffig) . ({k,-K)H 


— sine 


sine z — 


^ if |*| >0 

1 - if jzj -• 0 


(133) 


(134) 
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C Series Solutions of the Helmholtz Equation 
C.l Inner Products 

The inner product of two functions is defined as 

<u,v> = hJv u * vdI> (135) 

where T> indicates the domain of the functions. In one dimension, \D\ is the length of the 
interval; and, in two dimensions, it is the area. Because the first variable is conjugated, the 
average of a product over space is 

uv = (u*,v) (136) 

The average of the magnitude squared of a complex function is 

M 2 = (u,u) (137) 

and the average of the square of a function is 

u 2 = (u*,u) (138) 

The average of a squared function is complex, but the average of the magnitude squared is 
real. 

The inner product of two modes is defined as 

1 r+W/2 r+H/ 2 

(4>mn(y,z),<f> m ' n '(y,z)) = — — / / 4>* mn {y,z)4> m . n i(y,z)dydz (139) 

n vv J-W/2 J-H/2 

where the domain is HW . Because the two-dimensional modes are products of one-dimensional 
modes, 

{$mni fim'n ' ) = (XmiXm) (140) 

If we assume discrete sets of eigenvalues k ym and k zn , then the corresponding modal 
functions Xm(y) and i> n (z) are not orthogonal with this inner product; however, the complex 
conjugates of the modes form a reciprocal basis. Inner products of the conjugate modes and 
the modes are then proportional to Kronecker delta functions: 


(X*m{y),Xm>(y)) = Xm^rnm' 

(141) 

(^n( 2 )^n'(z)) = $ 2 S nn r 

(142) 

(^mniy i z), <f>m'n'{yt z)) — Xm^’n^mm'^nn' 

(143) 


where the overbar symbol means the average over the spatial dimensions of the function. 

C.2 Point Source in the Interior 

The three-dimensional wave equation with a point source of volumetric strength S is 

V 2 p s + k 2 p s = LopSS(x - x 0 )S(y - y 0 )S(z — z 0 ) (144) 

The pressure p s can be expressed in a modal series, as in equation (11) in terms of the basis 
functions 4> mn (y,z), which satisfy the wall boundary conditions: 

CO 

p a (x,y, Z ) = pc 2 53 Gmn<f>mn(y,z)e ,kxmnlX ~ X o| (145) 
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The combination of equations (144) and (145) gives 


-pc 2 E *) ( S + * 2 rn » ) e ’ 1 '""' 1 "" 1 

m,n=0 ' 


= u:pSti(x - x 0 )6(y - y 0 )6(z - xo), x ± xo 


(146) 


If we utilize the reciprocal basis function to form inner products and integrate over the 
6(x — x o) function, we obtain an equation for the coefficients Q mn • This leads to the condition 



( s ) 

Xm(y o) 

0n(2 O ) 

\HWcJ 

L Am J 



(147) 


C.3 Point Source in the Left Boundary 

Let the left boundary be a hard wall with a point source of strength S w at a point (y 0 ,z 0 ). 
Then, the admittance in the general boundary equation (2) is zero, and 


u 3 = S w S(y - y 0 )6(z - z 0 ) 
The velocity is expanded in a modal series as 


(148) 


h,{x,y,z) = C E x)e' kxmn * 

m.n—O 

By evaluating this expression at x — 0,we obtain the coefficients of the expansion 


Ut t 


(149) 


/ s w \ 

'Xm(yo) 

0n(2o) 

\HWc) 

L Am J 
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The pressure coefficients are given by the modal impedances 


V S+ = Zrr,„ 

r mn *-"mn 


/ \ 

Xm(yo) 

$n{z o) 

\HWc) 

V 2 

L Am J 

. $1 . 


(151) 


These coefficients differ from those for a source in the interior only by the factor —if 2. Note 
that a volumetric source in the wall is twice as effective as a volumetric source in the interior. 


C.4 Transmitted Power and Intensities for a Point Source 

The acoustic intensity at a point (x, y , z) is given by 

I x (x. y , z) = 29fc[pt£ m ] (152) 

where the superscript asterisk denotes the complex conjugate and 9? is the real part. The 
power transmitted in the axial direction is given by 

f H/2 f W/2 

WJx) = / / I(x,y, z)dydz = 2HW (u,p) (153) 

J-H/2 J-W/2 
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The inner product (u,p) is 


{u,p) = pc 3 

m,n=0 m',n'=0 


’mnM m t n re 


*( ^xmn / /) J^ — ^ol 


(Xrn'i Xrn) (0n'i0n) 


In the case of hard walls, the modes are orthogonal, and the inner products of modes 
Kronecker delta functions. The power formula then simplifies to a double sum: 


OO 

(' u,p)=pc 3 Vmn^ n e~ 2 ^ kimn) {x ~ x ° l 

771,71 = 0 


( 154 ) 

become 

( 155 ) 
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